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We study the zero-temperature ground-state (gs) phase diagram of the spin-^ anisotropic pla- 
nar pyrochlore (or crossed chain) model using the coupled cluster method (CCM). The model is 
equivalently described as a frustrated Ji~-J2 antiferromagnet on the two-dimensional checkerboard 
lattice, with nearest-neighbor exchange bonds of strength Ji > and next-nearest-neighbor bonds 
of strength J2 = kJi > 0. Using various antiferromagnetic (AFM) classical ground states as CCM 
model or reference states we present results for the gs energy, average on-site magnetization, and 
the susceptibilities of these states against the formation of plaquette valence-bond crystal (PVBC) 
and crossed-dimer valence-bond crystal (CDVBC) ordering. Our calculations show that the AFM 
quasiclassical state with Neel ordering is the gs phase for n < « 0.80 ± 0.01, but that none 
of the fourfold set of AFM states that are selected by quantum fluctuations at 0{l/s) in a large-s 
analysis (where s is the spin quantum number) from the infinitely degenerate set of AFM states that 
form the gs phase for the classical version (s — ^ 00) of the model (for k > 1) survives the quantum 
fluctuations to form a stable magnetically-ordered gs phase for the s = | case. We show that the 
quasiclassical Neel state becomes infinitely susceptible to PVBC ordering at or very near to n = , 
and that the quasiclassical fourfold AFM states become infinitely susceptible to PVBC ordering at 
K = ~ 1-22 ± 0.02. In turn, we find that these states become infinitely susceptible to CDVBC 
ordering for all values of k above a certain critical value at or very near to k = . Our calculations 
thus indicate a Neel-ordered gs phase for k < , a PVBC-ordered phase for < k < , and 
a CDVBC-ordered phase for n > Hc2- Both transitions are likely to be direct ones, although we 
cannot exclude very narrow coexistence regions confined to 0.79 < < 0.81 and 1.20 < k < 1.22 
respectively. 

PACS numbers: 75.10.Jm, 75.50.Ee, 75.40.-s, 75.10.Kt 



I. INTRODUCTION 

The theoretical study of two-dimensional (2D) frus- 
trated quantum antiferromagnets has been strongly mo- 
tivated by the fact that such quantum spin models of- 
ten describe well the properties of real magnetic mate- 
rials of great experimental interest. These models have 
also become of huge current interest because the inter- 
play between frustration and quantum fluctuations seen 
in them can produce, even at zero temperature (T = 0), a 
wide variety of fascinating quantum phases ranging from 
those with quasiclassical ordering to valence-bond solids 
and spin liquidsii"— They have thus become paradigms of 
systems that may be used to study quantum phase tran- 
sitions between quasiclassical phases showing magnetic 
order and magnetically disordered quantum phases. 

Some of the parameters that determine which type 



of ordering occurs include the lattice geometry, the di- 
mensionality D of the system, the spin quantum number 
s of the atoms situated on the lattice sites, the num- 
ber and range of the magnetic bonds, and the degree 
to which bond frustration of either the geometric or dy- 
namical kind is present. New impetus for the study of 2D 
quantum spin-lattice models comes from recent propos- 
als to realize them experimentally with ultracold atoms 
trapped in an optical latticei^ The particularly exciting 
scenario thus opens of being able to tune the competing 
bond strengths and thus to investigate experimentally 
the ensuing quantum phase transitions and their dynam- 
ics. 

One of the prime theoretical interests in frustrated 
quantum magnets lies in the possibility that they might 
exhibit quantum disordered states and/or spin-liquid be- 
havior. Among the most highly frustrated, and hence 
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most promising, candidate systems in this regard are 
those that are composed of tetrahedra coupled into two- 
dimensional (2D) or three-dimensional (3D) lattice net- 
works. Prominent among the latter are the pyrochlores, 
whose basic structure is one of vertex-sharing tetrahe- 
dra. Indeed, experiments on such s = ^ pyrochlores as 
Y2lr207 do seem to show evidence for a quantum spin- 
liquid state 1^ 

In order to reduce the complexity of the 3D pyrochlore 
lattice, but without diminishing the magnetic frustration, 
one may project the 3D vertex-sharing lattice of tetrahe- 
dra onto a 2D plane. Each tetrahedron comprises four 
spins at its vertices, with each of its six edges or links 
representing an interaction of the Heisenberg antiferro- 
magnet (HAFM) form. Each such tetrahedron is thus 
mapped to a square with spins at its vertices and with 
sides representing antiferromagnetic (AFM) bonds, but 
now with additional AFM links across its diagonals. Such 
a pattern is repeated in the vertex-sharing arrangement 
shown in the checkerboard pattern of Fig. [TJ Although 
this 2D projection of the 3D pyrochlore structure pre- 
serves its vertex-sharing structure, the symmetry in the 
3D structure between the six bonds on each tetrahedron 
is lost in the 2D projection since the two diagonal bonds 
of each crossed square are now inequivalent to the four 
bonds on the sides of the square. This subsequent reduc- 
tion in the symmetry is thus consistent with considering 
an anisotropic Heisenberg model on the 2D checkerboard 
lattice in which the AFM exchange interactions along the 
sides of the squares (with strength Ji > 0) are generally 
different in strength from those along the diagonals of the 
crossed squares (which have a strength J2 > 0), as shown 
in Fig. [TJ The resulting frustrated model is thus called 
the anisotropic planar pyrochlore. Alternative names are 
the anisotropic checkerboard HAFM, the J1-J2 checker- 
board model, and the crossed chain model. 

Although, the spin-i anisotropic planar pyrochlore has 
been studied by a large number of authors^"— the struc- 
ture of its full phase diagram still remains unsettled and 
contentious, especially for larger values of the frustra- 
tion parameter, k = J2/J1 ^ 1, as we discuss more 
fully below in Sec. [Ill Various methods have been ap- 
plied to the model for different regions of the parame- 
ter space for the variable k. These include semiclassical 
(s ^ 1) analyses,— large- expansions of the Sp(Af) 
model^ii^ii^ high-order cluster-based strong-coupling se- 
ries expansion (SE) techniques^ ' ^^1^^ using a continu- 
ous unitary transformation generated by the flow equa- 
tion method of Wegner,— a real-space renormalization 
techniquei^ using the contractor renormalization method 
of Morningstar and Weinstein^ an easy-axis general- 
ization of the 3D model,— a quasi-one-dimensional ap- 
proach (valid in the k ^ 1 limit) based on the ran- 
dom phase approximation backed up by a bosonization 
study,— techniques that combine renormalization group 
ideas with one-dimensional bosonization and current al- 
gebra methods,— exact diagonalization (ED) of small 
finite-lattice clusters j^iii^ii^ii^ a two-step density-matrix 



renormalization group method;^ and, very recently, a 
tensor network simulation^ based on infinite projected 
entangled pair states.— 

In this paper we use the coupled cluster method 
(CCM) of quantum many-body theory (see, e.g., 
Refs. |28l430l | and references cited therein) to study the 
spin-i J1-J2 Heisenberg model on the checkerboard lat- 
tice, in order to attempt to shed more light on it. The 
CCM has proven itself in many applications to frustrated 
magnetic systems to be capable of providing accurate es- 
timates of the quantum critical points marking the phase 
transitions between states of widely differing order (see 
e.g., Refs. [sot - lil p. In view of the continuing interest 
in the model and the controversy over its T = phase 
structure, especially at large frustration (k > 1), it seems 
appropriate and timely to bring to bear on the problem 
the proven power of the CCM. Since, as we shall see, we 
arc able to calculate to high orders in the relevant CCM 
approximation scheme, as discussed below in Sec. lIIIl we 
are able to present accurate results from a method of 
well-proven ability to deal with such strongly correlated 
and highly frustrated systems. 

We now briefly outline the structure of the remainder 
of the paper. In Sec. HI] the model itself is first described 
and discussed. The CCM formalism is then briefly out- 
lined in Sec. IIIII before we present and discuss our CCM 
results in Sec. |lVl Finally, we conclude in Sec. |V] with a 
summary of our findings and a comparison of our results 
for the model with those from other methods. 



II. THE MODEL 

The Hamiltonian for the anisotropic checkerboard- 
lattice model considered here is given by 

iJ = Jl ^ Sj • Sj + J2 ^ Si ■ Sfe . (1) 

where the index i runs over all sites of a square lattice, 
index j runs over all nearest-neighbor (NN) sites to site 
i, and index k runs over all next- nearest- neighbor (NNN) 
sites to site i on a checkerboard pattern such that alter- 
nate square plaquettes have either two NNN (diagonal) 
bonds or none, as shown in Fig. [1] The sums over 
and {{i,k)) count each pairwise bond once and once only. 
Each site i of the lattice carries a particle with spin s = ^ 
and a spin operator Sj = (sf , , sf ). 

The lattice and exchange bonds of the anisotropic 
checkerboard-lattice model are shown in Fig. [T] We may 
alternatively view the model as comprising crossed (diag- 
onal) sets of chains on which the intrachain exchange cou- 
pling constant is J2 , coupled by (vertical and horizontal) 
interchain exchange bonds of strength Ji. We assume 
here that both bonds are antiferromagnetic (AFM) in 
nature (i.e., have positive exchange coupling constants) 
and hence frustrate one another. The model thus inter- 
polates continuously between the isotropic HAFM on the 



3 

















































< 1 















1'' 


































< 



A A 



Y V 



(a) Neel 



(b) striped 



(c) Neel* 



FIG. 1: (Color online) The Jx-J^ checkerboard model (with J\ = 1), showing (a) the Neel state, (b) the (columnar) 
striped state and (c) one of the two Neel* states. The NN J\ bonds are shown as solid (black) lines and the NNN J2 
bonds arc shown as dashed (blue) lines. The arrows represent the orientations of the spins on each lattice site for 

each of the three states shown. 



square lattice (when k = J2/«/i = 0) and decoupled one- 
dimensional (ID) isotropic HAFM chains (when k — >■ 00). 
In between, at k = 1, we have the isotropic HAFM on 
the checkerboard lattice that is a 2D analog of the 3D 
isotropic pyrochlore HAFM. Henceforth, without loss of 
generality, we set Ji = 1 in order to set the energy scale. 

The classical ground-state (gs) phase for this model for 
K < 1 is the Neel state shown in Fig. [Ua), in which every 
column and row exhibits Neel AFM ordering, • • • titi 
• • • , and consequently the ordering along each diagonal 
is ferromagnetic (FM), i.e., where all the spins are aligned 
parallel to one another. The Neel state has an energy per 
spin given by E'^^/N = s^(— 2Ji + J2). For k > 1 there is 
an infinitely degenerate family of collinear gs phases in 
which every diagonal exhibits Neel AFM ordering, but 
where every diagonal, each of which is connected by Ji 
bonds to two other crossed diagonals, can be arbitrarily 
moved along its own direction. These states all have the 
same energy per spin of E'^^/N = — J2, independent of 
the exchange coupling Ji . The classical phase transition 
is clearly at KcI = 1 (</i > 0). 

Among this infinitely degenerate family of classical 
states for k > 1 are the so-called (columnar) striped 
state shown in Fig. [IJb) and the Neel* state shown in 
Fig. [TJc). The columnar (row) striped states have FM 
ordering along columns (rows) but AFM Neel ordering 
along rows (columns) . The Neel* state has doubled AFM 
ordering, • • • tt4-itt4-i ' • • : along every row and column. 
Thus, the single-site spin t or J, of the usual Neel state is 
replaced in the Neel* state by the two-site unit tt or ii- 
Like the striped state, the Neel* state is also doubly de- 
generate (for a given direction of the Neel vector) , since 
the roles of the rows and columns may be interchanged in 
Fig.HJc) (or, equivalently, the two-site unit of up or down 
spins may be chosen along rows as well as columns). 

Compared to the classical (s — ?> 00) version of the 
anisotropic checkerboard model, the s = | case is really 
only well established at the three points k = 0, k = 1, and 
K — !> 00. For the square-lattice HAFM (k = 0) almost all 
methods concur that the classical Nccl AFM long-range 
order (LRO) is not destroyed, although the staggered 



magnetization is reduced from the classical value of 0.5, 
and the excitations are gapless, integer-spin magnons. By 
continuity it is expected that the Neel order will persist 
as the frustrating J2-bonds are turned on, out to some 
critical value , at which the Neel staggered magneti- 
zation goes to zero. 

There is also a broad general consensus from a vari- 
ety of methods that at the isotropic point {k ~ 1) the 
gs phase of the s = ^ checkerboard-lattice HAFM is a 
plaquette valence-bond crystal (PVBC) with quadrumer 
LRO on isolated spin-singlet square plaquettes, and 
with gapped integer-spin excitations (that are confined 
spinous). It is still an open question as to whether there 
is a direct (first-order in the Landau-Ginzburg scenario) 
transition at k = between the states with Neel and 
PVBC order, or whether there is an intermediate coex- 
istence phase with two different order parameters. Such 
a phase could have continuous Landau-Ginzburg transi- 
tions to both the Neel and PVBC phases. The possibility 
of such coexistence regions occurring between Neel and 
valence-bond solids has been discussed in great detail 
both in a general context in Ref. [IH for various spin- 
lattice models, and in Ref. [2^ in the specific context of 
the present model. Again, by continuity, we expect that 
the PVBC order will persist to values of k out to some 
critical value > 1, at which point the PVBC order 
vanishes. 

Lastly, at the k — > cx) limit of the s = ^ anisotropic 
checkerboard model we have the well-known and ex- 
actly soluble case of decoupled ID HAFM chains. Such 
ID spin-i chains have a Luttinger spin-liquid gs phase, 
with a gapless excitation spectrum of dcconfincd spin-i 
spinous. 

The most unsettled part of the phase diagram for this 
model is the region n > Kc^, where various predictions 
have been given. For example, it has been arguedii that 
the ID Luttinger behavior of the k ^ 00 limit might 
be robust against the turning on of interchain (Ji) cou- 
plings, so that the chains continue to act as decoupled. 
Such a 2D spin-liquid gs phase provides an example of 
a so-called sliding Luttinger liquid4^"— Numerical evi- 
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dence for such a spin-liquid phase at large values of k 
in the present model was also found from ED studies on 
samples of up to iV = 36 spins.— 

Alternatively, by making a more careful analysis of the 
relevant terms near the ID Luttinger liquid fixed point, 
it was shown later— that the original predictionii of a 
sliding Luttinger liquid was wrong, and the same authors 
suggested that the correct gs phase in the large- k limit 
is the so-called gapped crossed dimer phase, where the 
system spontaneously dimerizes with a staggered order- 
ing of dimers along the J2 chains (i.e., along the diag- 
onals in Fig. [1]). Support for the crossed-dimer phase 
has come from series-expansioi*^^ and two-step density- 
matrix renormalization group method studiesi^ We dis- 
cuss this phase further in Sec. IIVI below. 

Finally, one may wonder whether any of the infinitely- 
degenerate set of classical (s — > 00) ground states for 
K > 1 may survive the quantum fluctuations present in 
the s = 5 model, and, if so, whether the classical de- 
generacy may be lifted by the well-known order by dis- 
order nicchanismi^ A semiclassical (s ^ 1) analysisii^ 
has shown that quantum spin-fluctuations induce a LRO 
that breaks the fourfold rotational symmetry of the lat- 
tice, and that to 0(l/s) the fourfold degenerate set of 
states comprising the striped state of Fig. [Tfb) and the 
Neel* state of Fig.[ljc) (plus their two counterparts where 
rows and columns arc interchanged) become energetically 
favored as the gs phase over the remainder of the infinite 
classical set. A very recent tensor network simulatioii^l 
of the spin-i model finds that, contrary to essentially all 
other calculations on this model, this fourfold degenerate 
state survives to be the quantum gs phase for all values 
of the frustration parameter above that at which PVBC 
order disappears (k > Kca). These authors also argue 
that, although their numerical program is unable to dis- 
tinguish between the energies of the striped and Neel* 
states in the quantum (s = ^) model, the striped phase 
will emerge as the actual gs phase in practice because of 
its greater robustness against small perturbations to the 
Hamiltonian. 

Other analyses2£ have, however, shown that, the Neel* 
state might, in one possible scenario, intervene as an in- 
termediate gs phase between the two (i.e., the plaquette 
and crossed-dimer) valence-bond solid phases. In such 
a scenario the transition between the PVBC and Neel* 
phases is shown to be able to proceed via a continuous 
0(3) transition, while that between the crossed dimer 
and Neel* phases will be either a direct (first-order in 
the Landau-Ginzburg scenario) one or will proceed via 
an intermediate coexistence phase showing both types of 
ordering (i.e., both Neel* spin ordering and crossed dimer 
bond modulation). 

In view of the considerable lack of agreement about 
the gs phase diagram for the spin-i J1-J2 model on the 
checkerboard lattice we now present results for it in the 
present paper from high-order CCM calculations. 



III. THE COUPLED CLUSTER METHOD 

The CCM (see, e.g., Refs. [H-ll^ and references cited 
therein) is one of the most powerful and universally ap- 
plicable quantum many-body techniques. It has been 
applied successfully to many quantum spin-systems (see 
e.g., Refs. (30l - l4^ and references cited therein). The 
method is particularly suitable for investigating highly 
frustrated quantum magnets, for which other alterna- 
tive methods may be of limited usefulness. For in- 
stance, quantum Monte Carlo (QMC) techniques are of- 
ten severely restricted by the well-known "minus-sign 
problem," which is ubiquitous in highly frustrated quan- 
tum magnets. On the other hand, ED methods are often 
too restricted by the relatively small size of the largest 
lattices that can be handled with given computational re- 
sources to be able to sample accurately the often subtle 
ordering present. 

Wc briefly describe the CCM formalism here and we 
refer the interested reader to the literature (and see, 
e.g., Refs. [30l - [39| and references cited therein) for fur- 
ther details. The implementation of the CCM always 
begins with the choice of a suitable reference or model 
state. It is usual, but by no means vital, to choose a 
classical gs phase as the model state |<I>). Hence, for 
the present anisotropic checkerboard model, we choose 
the Neel state, the striped state and the Neel* state as 
our CCM model states. From the discussion in Sec. [Tl] 
above we expect that the Neel state is likely to provide 
a good candidate CCM model state in the region k < 1, 
while the striped and the Neel* states are expected to 
be suitable candidates for k > 1. We choose only the 
latter states out of the infinitely degenerate set of classi- 
cal states in the k > 1 regime since, as discussed above, 
this fourfold set of states is selected by the order by dis- 
order mechanism at the 0(l/s) level in a quasiclassical 
expansion in powers of l/sji^ at which order they remain 
degenerate in energy. 

The CCM then incorporates the multi-particle correla- 
tions present in the exact quantum gs phase under inves- 
tigation on top of the chosen model state in a systematic 
hierarchy of approximations for the correlation operators 
S and 5' which parametrize the gs ket and bra wave func- 
tions as 

|«') = e^|$); (^1 = ($|S'e-'^. (2) 
The correlation operators are written as 

S = Y,SiC+; S = Y,SiCY; V/ ^ 0, (3) 

where Cq = 1, the identity operator, / is a set-index 
describing a set of single-particle configurations, and 
and = (C/)^7 foi' 7^ 0, are Hermitian-conjugate 
pairs of multi-particle creation and destruction operators 
defined with respect to the model state 1$) considered as 
a generalized vacuum state. They are thus required to 
satisfy the conditions ($|C/ = = C7|$);V/ ^ 0. They 
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form a complete set of mutually commuting many-body 
creation operators in the Hilbert space, defined with re- 
spect to |4>) as a cyclic vector. The states are normalized 
such that (^-l^-) = (^l^-) = ($1$) = 1. 

For spin-lattice systems it is convenient to choose a 
set of local coordinate frames in spin space such that on 
each lattice site the spin in each model state points in 
the downward (negative z direction). Such rotations ob- 
viously do not affect the basic SU(2) spin commutation 
relations, but they have the simplifying effect that the op- 
erators are transformed into multi-spin raising opera- 
tors that can be expressed as products of single-spin rais- 
ing operators, C/ = s^s^ • • • s^, where Sj' = Sj+is^. 

The gs energy is evaluated in terms of the correlation 
coefficients {Si} as E = = {<P\e-''^ He^\<P)] and 

the average on-site magnetization AI in the rotated spin 
coordinates is evaluated cquivalently in terms of the co- 
efficients {Si, Si} as AI = -^(^| E^i s|l*>- Thus, M 
is simply the usual magnetic order parameter. 

The complete set of unknown ket- and bra-state cor- 
relation coefficients {<S/,5/} is evaluated by setting the 
energy expectation value H = {'^\H\'^) to be a mini- 
mum with respect to all parameters {Si,Si;yi ^ 0}. 
This produces the coupled set of nonlinear equations for 
the ket-state (creation) correlation coefficients {Si} via 
($|C7c-'5He'^|$) = 0;V/ ^ 0; plus the coupled set of 
linear equations, {^S{(r^Ue^ - E)Cl\<P) = 0;V/ 7^ 0, 
which are used to solve for the bra-state (destruction) 
correlation coefficients {Si}. 

If it were possible to consider all creation and anni- 
hilation operators Cf and respectively, i.e., all sets 
(configurations) of lattice sites, in the CCM correlation 
operators S and S respectively, one would in principle ob- 
tain the exact eigenstatc of the system belonging to any 
symmetries imposed by the model state (and the config- 
urations that are perhaps also accordingly selected))^ Of 
course, however, it is necessary in practice to use approx- 
imations schemes to truncate the expansions of S and S 
in Eq. ([3]). In that case the approximate results for the 
gs energy E and the magnetization M will depend on the 
choice of model state. 

For the case of s = i systems, as considered here, we 
normally use the well-tested localized LSUBto truncation 
scheme which takes in at the mth level of approximation 
all multi-spin correlations in the CCM correlation oper- 
ators over all configured regions on the lattice defined by 
m or fewer contiguous sites. A configuration of m sites 
is considered to be contiguous if every site in the con- 
figuration is adjacent (in the NN sense) to at least one 
other site in the configuration. Clearly, as m 00, the 
LSUBm approximation becomes exact. For the present 
checkerboard model, we use the CCM and the LSUBm 
scheme with m < 10 for the three model states shown in 
Fig.[TJ For the LSUBm configurations we assume the fun- 
damental checkerboard geometry to define the LSUBm 
sequences, and hence treat both the pairs of sites con- 
nected by Ji bonds and those connected by J2 bonds as 



TABLE I: Number of fundamental configurations, Nf, 
for the checkerboard geometry for the spin-i J1-J2 
checkerboard model (Ji = 1), using the Neel, striped, 
and Neel* states as CCM model states. 



Method 


Nf 


Neel 


striped 


Neel* 


LSUB4 


27 


54 


79 


LSUB6 


632 


1225 


2441 


LSUB8 


21317 


41324 


86590 


LSUBIO 


825851 


1598675 


3373495 



being contiguous sites. Table U shows the number Nf 
of such distinct (i.e., under the symmetries of the lattice 
and the model state) fundamental spin configurations for 
each of the three model states that we use for our spin-i 
J1-J2 checkerboard model. 

It is clear that Nf rises rapidly with the truncation 
index m. For example, for the Neel* state that is used in 
this study as one of the CCM model states, the LSUBIO 
approximation contains 3373495 distinct spin configura- 
tions. This is the highest LSUBm level that we can reach 
here using the Neel* state as our model state, even with 
massive parallelization and the use of supercomputing re- 
sources. It takes us approximately 1 h computing time 
using massively parallel computing with 3000 processors 
simultaneously to solve the corresponding coupled sets 
of CCM bra- and ket-state equations, to obtain a single 
data point for a given value of J2, with Ji = 1^ 

We note that if, instead of using the checkerboard ge- 
ometry, we were to use the square-lattice geometry (i.e., 
with NN pairs defined only by Ji bonds), the number of 
fundamental configurations Nf would obviously be fewer, 
at the same level m, than in the checkerboard geometry. 
In turn this could perhaps enable us go to higher LSUBm 
orders for given computational power. However, this ad- 
vantage is completely outweighed by the disadvantage 
that the LSUBm sequences for both E/N and A4 then 
show a marked staggering behavior in m = 2k, depend- 
ing on whether k is even or odd. This is clearly due to the 
fact that the full LSUBm sequence does not then prop- 
erly reflect the checkerboard symmetries. It is quite sim- 
ilar to the odd and even staggering behavior in index m 
for LSUBm approximations on simple (dynamically un- 
frustrated) models, which has been reported elsewherei^ 
Any such staggering effect makes extrapolations (for the 
full sequence) of the sort we now discuss more compli- 
cated and less robust. 

Thus, as a final step we need to extrapolate the raw 
CCM data from our LSUBm approximations to the ex- 
act (m — >■ 00) limit. In the absence of any staggering 
effects of the sort described above, we use the well-tested 
extrapolation scheme 

E{m)/N — Go + flim^^ + a2m~^ , (4) 

for the gs energyj^^— For the magnetic order param- 
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FIG. 2: (a) The extrapolated CCM LSUBcx) results for the gs energy per spin, e = E/N, versus J2 for the striped 
and Neel* phases of the spin-i J1-J2 Heisenbcrg antiferromagnet on the checkerboard lattice (with Ji = 1), using 
the LSUBm resuhs with m = {4, 6, 8, 10} and Eq. (g]); (b) the energy difference, Ae = e^*"?'''^ - e^°'='* versus J2 of 
the two phases shown in (a) using LSUBm approximations with m = {4, 6, 8, 10} and also using the corresponding 
separate LSUBoo results for both phases from Eq. (|4]) using m = {4, 6, 8, 10}. 



eter, M, we use the schemes 

M{m) = bo + him~^ + 62™"^ : 
for non-frustrated spin systems)^"— and 



M{m) = Co + bim + a2m 



,-3/2 



(5) 



(6) 



for highly frustrated spin systems.— i^>^2ii^ We have 
performed separate extrapolations using data sets with 
TO = {4,6,8,10}, TO = {6,8,10}, TO = {2,4,6,8}, and 
m ~ {4,6,8}. They yield very similar results in each 
of the cases reported below, which gives credence to our 
results and demonstrates their robustness. 
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IV. RESULTS AND DISCUSSION 

We now present our CCM results for the spin-^ 
anisotropic checkerboard model, using each of the three 
states shown in Fig. [1] as model states. We first show 
in Fig. [2]^a) the extrapolated LSUBoo gs energies per 
spin, E/N, of the phases obtained using the striped and 
Neel* model states. We recall that in the classical limit 
(s — 00) these two phases are degenerate and are the 
gs phase only for k > 1. Results are shown in Fig. ^a.) 
down to the lowest terminating values of k in each case for 
which real solutions exist for all of the LSUBto approx- 
imations used. It has been shown previousl y'^ '^i^^ that 
such termination points are a strong indication of the cor- 
responding quantum phase transition points that occur 
in the system under study. Figure [5Jb) shows the differ- 
ence in the energies of the two states in the approximate 
region where both CCM solutions exist. We see clearly 
that although the energy difference is small, the classical 
degeneracy is removed in favor of the Neel* state over 



FIG. 3: CCM results for the gs energy, E/N, for the 
Neel and Neel* states for the spin-i J1-J2 Heisenberg 
antiferromagnet on the checkerboard lattice (with 
Ji = 1) versus J2. The LSUBto approximations for 
TO = {4, 6, 8, 10} are shown together with the 
corresponding LSUBoo extrapolation from using Eq. 
with TO = {4,6,8,10}. 



the striped state for all values of k for which real CCM 
solutions for both phases exist. Nevertheless, bearing in 
mind the smallncss of the energy difference, we present 
some results below using both states as model states. 

In Fig. [3] we show both the LSUBto and the extrap- 
olated LSUBoo results for the gs energy, E/N, of both 
the Neel and Neel* phases. As noted briefly above, we 
now observe more clearly that each of the LSUBto en- 
ergy curves based on a particular model state termi- 
nates at some critical value of k (that itself depends on 
the LSUBto approximation used), beyond which no real 
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(a) Neel and Neel* states (b) Neel and striped states 

FIG. 4: CCM results for the gs magnetic order parameter, M, for the Neel and Neel* states for the spin-^ J1-J2 
Hcisenberg antifcrromagnet on the checkerboard lattice (with Ji = 1) versus J2; (a) using the Neel and Neel* states 
as model states, and (b) using the Neel and striped states as model states. In both cases the LSUBm 
approximations for m = {4, 6, 8, 10} are shown together with the corresponding LSUBoo(l) and LSUBcxd(2) 
extrapolation from using Eqs. ([5]) and © respectively, with m = {4, 6, 8, 10}. 



CCM solution can be found. Since the CCM LSUBm 
solutions require increasingly more computational power 
to obtain (to a given level of numerical accuracy) as the 
termination points are approached, it is computationally 
costly to determine the actual termination points to a 
high degree of accuracy. We note that in Fig. [3] results 
are shown for each LSUBm case down to values of k be- 
low which for the Neel* phase, and up to values of k above 
which for the Neel phase, real solutions based on the re- 
spective model state cease to exist. As noted above, in 
all cases the corresponding termination point at a given 
LSUBm level shown in Fig.[2]for the striped state is lower 
than that for the equivalent Neel* model state case. 

We note however that, as is usually the case, the CCM 
LSUBm results for finite m values for both the Neel and 
Neel* phases shown in Fig. [3] extend beyond the cor- 
responding LSUBoo transition points. For large values 
of m the LSUBm transition points are quite close to 
the actual quantum critical points (QCPs) where that 
phase ends. For example, the LSUBIO termination points 
shown in Fig.[3]are at « 0.88 for the Neel state and 
^NoGi- ^ 2 for the Neel* state. The CCM results show 
a clear intermediate regime in which neither of the qua- 
siclassical AFM states (Neel and Neel*) is stable. 

We now discuss the magnetic order parameter (viz., 
the average on-site magnetization), M , in order to inves- 
tigate the stability of the quasiclassical magnetic LRO. 
Our CCM results for M are shown in Fig. |4] for each of 
the Neel, Neel*, and striped phases. Our extrapolated 
results for M in the Neel phase are seen to be somewhat 
sensitive to whether we use the scheme of Eq. (O or that 
of Eq. ([6]). As we have indicated previously the scheme of 
Eq. (O is appropriate only for small values of J2 ■ For ex- 
ample, in the square-lattice limit, J2 = 0, we obtain the 
extrapolated LSUBoo(l) result M « 0.3069 from the use 



of Eq. © and the LSUBm values with m = {4, 6, 8, 10}. 
Very similar values are obtained with the alternative data 
sets m = {4, 6, 8} and m = {6, 8, 10}. We note that for 
the square-lattice HAFM no dynamic (or geometric) frus- 
tration exists and the Marshall- Peierls sign rule^ applies 
and may hence be used to circumvent the QMC "minus- 
sign problem." The QMC result^ M = 0.3070 ±0.0003, 
is thus extremely accurate for this limiting ( J2 = 0) case 
only. Our own CCM result using the scheme of Eq. ([5]) 
is thus in excellent agreement with it. By contrast, the 
extrapolation scheme of Eq. ([H]), which is appropriate for 
(highly) frustrated systems, gives a much poorer estimate 
of M w 0.275. 

The magnetization results show clear evidence for the 
melting of Neel order at a value k = < = 1, 
with results for that are very close to the corre- 
sponding termination point nf^'^^ discussed above. At 
such values k 1, where the system is highly frustrated, 
the extrapolation scheme of Eq. ^ is more appropri- 
ate, as we have indicated previously, and its use with 
the LSUBm data set m = {4, 6, 8, 10} gives us our first 
estimate for the quantum critical point (QCP) at which 
Neel order vanishes, « 0.796. Very similar results 
are found by using the alternative data sets m = {4, 6, 8} 
and m = {6,8,10}. Combining all these results gives 
the estimate ~ 0.80 ± 0.01. (By contrast, the use of 
the scheme of Eq. which is inappropriate in this frus- 
trated region near a QCP, gives a value Kci ~ 0.87±0.01.) 

The results in Fig. (Ha) for M using the Neel* state 
as CCM model state are seen to be qualitatively very 
different from those using the Neel state as model state. 
Indeed, all of the evidence from Fig. lUJa) is that M is 
either zero or very close to zero over the entire range for 
which CCM extrapolated LSUBcxd solutions exist using 
the Neel* state as model state. The more appropriate 
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extrapolation scheme of Eq. ([6]) in this regime with k > 1 
gives either negative values for M or positive values very 
close to zero over the entire range shown in Fig. |4j while 
even the inappropriate scheme of Eq. ([5]) gives only a very 
small and almost constant value of M « 0.08 over the 
same range. Our results in the high-frustration regime 
(k 3> 1) using the appropriate scheme of Eq. ^ are 
given extra credence by the fact that wc see clearly from 
Fig- HI that M — > rather accurately in the large k limit, 
which is the exact result for this limit where the model 
reduces to unlinked ID spin-^ chains. It is clear that at 
best the existence of the Neel* phase in the spin-i case 
is extremely fragile from this evidence. More likely, it is 
not the stable gs phase for any value of k, based on the 
results for M. 

For this reason we have repeated the calculations for 
M, but now using the striped state as CCM model state, 
even though we found it to have a slightly higher energy 
for all values of k than that of the Nccl* state. Results 
are shown in Fig.|4|b). Figures |4lja) and|4|b) show very 
similar results for M for the striped state and the Neel* 
state, except very near their corresponding termination 
points. All of the evidence so far is that neither state is 
the stable gs phase for any value of n. Since the results 
for the order parameter M are so similar for the Neel* 
state and the striped state, and since the former has a 
slightly lower energy, we henceforth restrict ourselves for 
larger values of k to use of the Neel* state as CCM model 
state. 

It is reasonably well established from earlier numer- 
ical studies using EDi^ and strong-coupling expansion 
techniques^i^ii^ that the gs phase of the spin-i HAFM 
on the checkerboard lattice (i.e., our model at the 
isotropic point J2 ~ Ji) is a. plaquette valence-bond crys- 
tal (PVBC) with long-range quadrumer order. Further 
evidence for such a valence-bond solid built from discon- 
nected 4-spin singlets comes from a "fermionic" SU(n) 
generalization of the SU(2) group in the large-n limit 
There is broad agreement from all this work that the 
PVBC phase comprises singlet plaquettes on the squares 
in Fig. [T] without crossed links, as shown in Fig.[SJ Thus, 
in order to get more information on the phase that oc- 
curs after the melting of Neel order at k = we now 
investigate the possibility that it might be a PVBC state 
of the sort shown in Fig. [SJ 

To do so we first consider a generalized susceptibility 
Xf that describes the response of the system to a pertur- 
bation described by a "field" operator F. A field term 
F = 6 Of is thus added to the Hamiltonian of Eq. ([T]). 
The energy per site in a given state, E{S)/N = e{5), is 
then calculated for the perturbed Hamiltonian H + F, 
and the susceptibility of the system to the perturbation 
F is defined a,s xf = — id'^e{6))/ {dS'^)\g^^. An instabil- 
ity of the state against the perturbation F is signalled by 
a zero point of Xf^ or, equivalently, by a divergence of 
Xf- In our case we first use the CCM to calculate xf, 
using a specific model state, in various LSUBm approx- 
imations. Although rather less empirical experience is 
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FIG. 5: (Color online) Left: CCM results for the inverse 
plaquette susceptibility, 1/xp, versus J2, using the Neel 
and Neel* states as model states, for the spin-i Ji~J2 
Heisenberg antiferromagnet on the checkerboard lattice 
(with Ji = 1). The LSUBm approximations for 
m = {4, 6, 8, 10} are shown together with the 
corresponding LSUBoo(l) and LSUBc>o(2) results from 
using the extrapolations schemes of Eqs. ([7]) and © 
respectively, with m = {4,6,8,10}. Right: The 
perturbations (fields) F = 60p for the plaquette 
susceptibility Xp- Thick (red) and thin (black) lines 
correspond respectively to strengthened and weakened 
NN exchange couplings, where Op = J2{i j) ^ij^i ' ^j, 
and the sum runs over all NN bonds, with aij = +1 and 
— 1 for thick (red) and thin (black) lines respectively. 



available for the m — > cxo extrapolation of the CCM data 
for Xf than for other quantities such as the gs energy E 
or the order parameter M, we have found previously^ 
that the same extrapolation used for the gs energy [i.e., 
XF(m) = rfo + dim^^ + c?2W~^] fits the data most accu- 
rately, at least in regions not too close to a divergence 
of the susceptibility. We also saw previously^ that a 
corresponding extrapolation of the inverse susceptibility, 

Xp^ {m) = xo + xim^'^ + X2m^'^ , (7) 

gave very consistent results that agreed very closely with 
those from the corresponding above extrapolation of xf, 
except again in regions close to where Xf^ ~^ 0- Since, 
as we will see below, we will be especially interested in 
precisely such regions over a large range of values of k, 
we also use the fitting function, 

Xf\^) = yo + yim-y^ ■ (8) 

To calculate the susceptibility, Xp^ of our system 
against PVBC ordering we thus set the operator to 
Op as illustrated in the right panel and the caption of 
Fig. O The perturbation field F thus breaks the trans- 
lational symmetry of H. Wc show CCM results in the 
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left panel of Fig. [5] using both the Neel and Neel* states 
as model states. We first observe that for smaller val- 
ues of J2 (i.e., on the Neel side) the two extrapolations 
agree very closely, even near the point at which Xp^ goes 
to zero. Thus, the extrapolated inverse plaquette sus- 
ceptibility using the LSUBm data set m = {4,6,8,10} 
(m = {6, 8, 10}) vanishes on the Neel side at k w 0.843 
(k « 0.833) using the extrapolation scheme of Eq. (O, 
and at K « 0.820 {k w 0.775) using the extrapolation 
scheme of Eq. ([5]). Since the exponent 2/2 in Eq. ([5]) falls 
rather sharply to a value close to 1 near the point at 
which vanishes, the best estimate for this point is 
more likely to come from the extrapolation scheme of 
Eq. dS]) than from that of Eq. 

Combining all these results gives our best estimate 
of K « 0.79 ± 0.03 for the point at which the Neel 
phase becomes susceptible to PVBC ordering. This 
is in very good agreement with the above estimate of 
Kci ~ 0.80 ± 0.01 at which Neel LRO disappears as mea- 
sured by our results for the order parameter M. Thus, 
our results show no evidence at all for a coexistence re- 
gion in which Neel and PVBC ordering are both present, 
such as has been suggested might occur^^ although we 
cannot exclude the possibility of a very narrow region of 
coexistence confined to the region 0.79 ^ k < 0.81. Our 
findings are in agreement with ED results^^ for the same 
spin-i anisotropic planar pyrochlore model that reach 
the conclusion that, if present at all, any such coexistence 
region is very narrow indeed. As has been discussed in 
detail elsewhere;^ the QCP at between the Neel and 
PVBC phases is both forbidden as a continuous transi- 
tion within standard Landau-Ginzburg theory and docs 
not seem either to be a viable candidate for a deconfined 
(continuous) transition. The shape of the magnetization 
curves in Fig. S] on the Neel side, which show a rapid de- 
crease near , is perhaps more indicative of a first-order 
transition, as we have observed previously?^ although 
such evidence should not be regarded as conclusive. 

We also observe from Fig. [5] that with the Neel* state 
used as our CCM model state the two extrapolations of 
Eqs. ([7]) and ^ for the inverse plaquette susceptibility, 
Xp^i agree quite closely at larger values of k but diverge 
slightly at smaller values, where itself becomes small. 
As K — cx) we observe that the exponent y2 in Eq. ([5]) 
appears to approach the value 1.5 [rather than 2 as in 
Eq. ([7])], and then drop to a value close to 1 as ap- 
proaches zero. For these reasons again, we expect the 
extrapolation of Eq. ([5]) to be more exact, especially in 
regions where Xp^ becomes small. Thus, the extrapo- 
lated inverse plaquette susceptibility using the LSUBm 
data set m ~ {4, 6, 8, 10} (m ~ {6, 8, 10}) vanishes on 
the Neel* side at k « 1.238 (k « 1.216) using the extrap- 
olation scheme of Eq. (|S]). Although, as we have seen, 
the Neel* state does not appear to exist as as a stable gs 
phase (since its order parameter M seems to vanish for 
all values of k, nevertheless the results using it as a CCM 
model state provide a robust basis for the calculation of 
Xp, and give an estimate ~ 1.22 ± 0.02 for the upper 



QCP at which PVBC order disappears. 

However, we are now led to the question of what is 
the actual gs phase of the model for larger values of frus- 
tration. K > Kc2, beyond the upper QCP (at k — Kc^) 
at which the PVBC phase ceases to exist as a stable gs 
phase. In the first place it is clear that such a QCP must 
exist since in the limit k ^ 00 one has the physics of 
decoupled HAFM ID chains, which are known to exhibit 
Luttinger spin-liquid behavior, and are typified by a gap- 
less excitation spectrum and spin-spin correlations that 
decay algebraically with inter-spin separation distance. 
It was arguedii that for large values of k, where the ID 
chains are weakly coupled, the gs phase might be a 2D 
sliding Luttinger liquid phase (characterized by the ab- 
sence of LRO and with massless deconfined spinous as 
the elementary excitations) that joined smoothly to the 
K — > 00 limit. It was later shown,— by a more care- 
ful analysis of the relevant terms near the Luttinger liq- 
uid fixed point of the independent ID spin chain, that 
this finding was incorrect. Instead, using techniques that 
combine renormalization group ideas and ID bosoniza- 
tion and current algebra methods, it was shown that 
in this large-K regime the gs phase is of spontaneously 
dimerized type, with a staggered ordering of dimers along 
the parallel chains (viz., the diagonals in Fig. [T]). Such 
a crosscd-dimer valence-bond crystal (CDVBC) phase, 
with twofold spontaneous symmetry breaking and no 
magnetic order, was independently confirmed^^ by a se- 
ries expansion (SE) technique based on the flow equation 
method. 

The CDVBC phase is illustrated in the right panel of 
Fig. |6l Similarly to what was done above for the PVBC 
susceptibility, Xpi we can now calculate the susceptibility, 
Xd of our system against CDVBC ordering by setting the 
perturbation operator Op to the operator Od illustrated 
in the right panel and the caption of Fig. [6l Since in the 
large J2 limit the energy scales linearly with J2 (as may 
clearly also be observed from Figs. HJa) and |31 we show 
our CCM results in Fig. IH] for the scaled inverse dimer 
susceptibility, J2 /xd as a function of J2 ■ 

Interestingly, in this case, the extrapolation scheme of 
Eq. ^ docs not fit the LSUBm data points at all well, 
and consequently gives a very poor fit. The reason be- 
comes quite evident when the extrapolation scheme of 
Eq. ^ is used instead. It is then observed that that the 
scaling exponent 2/2 approaches the value 0.75 for large 
values of J2, and even for smaller values near the QCP 
at K = only rises slightly to values that approach 1. 
It is clear that the extrapolated value of J2/Xd is small, 
slightly negative, and almost constant for all values of k 
shown, when the more flexible extrapolation scheme of 
Eq. ([H]) is used. Our results are consistent with the inter- 
pretation that the inverse dimer susceptibility is zero for 
all values k shown in Fig. [SJ namely those where we have 
real solutions to the equations pertaining to all of the 
CCM LSUBm schemes used. The actual lower termina- 
tion point is not easy to determine accurately in this way, 
as we have already mentioned above. However, it is quite 
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FIG. 6: (Color online) Left: CCM results for the scaled 
inverse crossed dimer susceptibility, J2/xd, using the 
Neel* state as model state, for the spin-i J1-J2 
Heisenberg antiferromagnet on the checkerboard lattice 
(with Ji = 1) versus J2. The LSUBm approximations 

for m ^ {4, 6, 8, 10} are shown together with the 
corresponding LSUBoo(l) and LSUBc>o(2) results from 
using the extrapolations schemes of Eqs. ([7]) and ([5]) 
respectively, with m = {4, 6, 8, 10}. Right: The 
perturbations (fields) F — SOd for the dimer 
susceptibility Xd- Thick (red) dashed and thin (blue) 
dashed lines correspond respectively to strengthened 
and weakened NNN exchange couplings, where 
Od = J2{{ik)) ^ifc^i ' ^k, and the sum runs over the NNN 
diagonal bonds of the checkerboard lattice, with 
= +1 and —1 for thick (red) dashed and thin (blue) 
dashed lines respectively. 



clear from Fig. [S] that the individual LSUBm curves all 
have a minimum at a value k ~ 1.20 ± 0.05, and our re- 
sults are thus consistent with the interpretation that the 
inverse dimer susceptibility is zero for all values k > ■ 
In that scenario there is a QCP between the PVBC and 
CDVBC phases at k = Kc2, and the CDVBC phase then 
persists out to the k ^ 00 limit of unlinked ID chains, 
which is itself a singular point. 

Our results give us essentially no information, how- 
ever, on the nature of the transition at Kc2- Starykh et 
al^ have concluded that a continuous transition between 
the PVBC and CDVBC phases is prohibited within the 
standard Landau-Ginzburg scenario of phase transitions. 
They argue further that the most probable alternative is 
a direct first-order transition, although a separate pos- 
sibility is again the existence of an intermediate coexis- 
tence phase with both valence-bond orderings that can 
then have separate continuous transitions to both the 
PVBC and CDVBC phases. From our results we clearly 
favor the former scenario, although we cannot exclude 
the possibility of a very narrow such coexistence phase 
confined to the region 1.20 ^ k ^ 1.22. 



SUMMARY AND CONCLUSIONS 



To summarize, we have investigated the gs properties 
and (T = 0) phase diagram of the frustrated spin-i an- 
tiferromagnetic J1-J2 model (with J2 = kJi ; Ji > 0) on 
the 2D checkerboard lattice, using the CCM carried out 
to high orders. 

In common with most other calculations on the model 
we find that the gs phase is an AFM Neel-ordered state 



for 



K < K. 



Cl J 



at which point the staggered Neel magne- 



tization vanishes. Our best estimate for this lower QCP 
is ~ 0.80 ± 0.01. This is in reasonable agreement, 
but probably more accurate than, a corresponding esti- 



mate of 



0.75 from an ED studj™ on samples of 



N = 16, 32, 36 spins. Since we calculate that the Neel- 
ordered state becomes susceptible to PVBC ordering at 
K « 0.79 ± 0.03 our results point then to a direct transi- 
tion from the Neel-ordered gs phase to a PVBC phase 
at K = , although we cannot exclude entirely the 
small possibility of a vey narrow coexistence region. This 
finding is in good agreement with that found from ED 
studies We estimate that any such coexistence region 
of AFM Neel ordering and quadrumer plaquette ordering 
is confined to the parameter range 0.79 ^ k < 0.81. Our 
estimate for the lower boundary at which the quadrumer 
PVBC ordering vanishes is in excellent agreement with 
the corresponding value k sa 0.80±0.01 from a high order 
SE calculation^^ that starts from the limit of uncoupled 
quadrumers. A recent tensor network simulation of the 
mode]3^ gives a somewhat higher value of « 0.88 for 
the QCP from the Neel gs phase to the PVBC gs phase. 

From our CCM calculations we estimate that the 
quadrumer order of the PVBC state vanishes at a higher 
QCP at K = ~ 1.22 ± 0.02. This is in reasonable 
agreement, although somewhat higher than the corre- 
sponding estimates k w 1.095 ± 0.035 from a high-order 
SE calculation^ and k « 1.11 from a tensor network 
simulation)^ 

Although the striped and Neel* states, which are the 
fourfold-degenerate quasiclassical gs phases at 0(l/s) in 
an expansion in powers of 1/s for k > provide excel- 
lent model states for CCM calculations at larger values 
of K > 1 in the sense of providing well-converged LSUBm 
results, we find that they are not stable gs phases for the 
spin-i model for any value of k. This is in sharp dis- 
agreement with the finding from a recent tensor network 
simulation^! that the striped and/or Neel* states form 
the stable gs phase for all values k > ■ By contrast we 
find from our CCM calculations that the Neel* state is 
susceptible to the formation of a CDVBC phase (with an 
inverse susceptibility that is essentially zero) for all values 
K > 1.20 ± 0.05. We conclude that the CDVBC state is 
thus likely to be the stable gs phase for k > , although 
we cannot exclude the possibility of a very narrow coex- 
istence regime confined to the range 1.20 ^ k ^ 1.22. 
We find no evidence at all for any region separating the 
PVBC and CDVBC phases where the Neel* state would 
form a stable gs phase, thereby ruling out one of the two 
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scenarios postulated by Starykh et al.^ although now 
in complete agreement with their other scenario for the 
gs phase diagram, on which we now provide accurate nu- 
merical results for the two QCP's at Kci and Kc2- 

Finally we note that the CCM has been used here with 
the model or reference states chosen as classical states 
built from independent-spin product states. For larger 
values of the frustration parameter, k > , we remark 
that these model states are able (with appropriate choice 
of additional perturbative terms in the Hamiltonian), to 
describe the susceptibilities of the system to form the true 
gs phasaes perfectly well, even though we have shown 
that the respective extrapolated order parameters with 
respect to these model states are essentially zero over 
the entire range. Nevertheless, it would be worthwhile 
to repeat the CCM calculations directly with dimer and 



plaquette valence-bond states. Indeed, just such a gen- 
eral CCM approach^ has recently been described, which 
combines exact solutions for dimer or plaquette valence- 
bond solid ground states with the computational imple- 
mentation described here^ that is based on independent- 
spin product model states. It would be interesting to use 
this formalism for the present model in order to confirm 
our results. 
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